Majorana fermions in vortex lattices 
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We consider Majorana fermions tunneling among an array of vortices in a 2D chiral p-wave su- 
perconductor or equivalent material. We calculate the amplitude for Majorana fermions tunneling 
between a pair of vortices and find that it depends on the background superconducting phase profile; 
it is proportional to the sine of half the difference between the phases at the two vortices. Using 
this result we study tight-binding models of Majorana fermions in vortices arranged in triangular 
or square lattices. In both cases we find that the aforementioned phase-tunneling relationship leads 
to the creation of super lattices where the Majorana fermions form macroscopically degenerate lo- 
calizable flat bands at zero energy, in addition to other dispersive bands. This finding suggests that 
tunneling processes in these vortex arrays do not change the energies of a finite fraction of Majo- 
rana fermions, contrary to previous expectation. The presence of flat Majorana bands, and hence 
less-than-expected decoherence in these vortex arrays, bodes well for the prospects of topological 
quantum computation with large numbers of Majorana states. 

PACS numbers: 03.67.Lx, 03.65.Vf, 71.10.Pm, 73.43.Jn 



The last two decades have seen a flurry of activity in 
the fields of topological states of matter and topological 
quantum computation(TQC)[T]. One candidate particle 
for TQC with the requisite exotic 'anyonic' characteris- 
tics is the Majorana fermion state in the vortex of a chi- 
ral p-wave superconductor [2j [3]. Such Majorana states 
are predicted to arise in the v = 5/2 fractional quantum 
Hall state [4], Sr2Ru04[5 J, topological insulator surface 
state [6] and in nanowires with spin-orbit coupling [7], to 
name but a few widely-reported possibilities. In 2D sys- 
tems there have been proposals to create multiple vor- 
tices in the host superconductor and to perform TQC 
by 'braiding' them, i.e, by physically moving the vortices 
around each other [2]. Practical implementations of TQC 
would require bringing together a large number of Majo- 
rana fermions and much attention has been paid to the 
splitting of their energy levels due to inter- vortex tun- 
neling. These level splittings would introduce decoher- 
ence on timescales inversely proportional to their mag- 
nitudes and place duration constraints on computation 
processes. In 2D systems, prospects for realizing large 
vortex arrays are bright [8 , potentially enabling the cre- 
ation of vortex lattices with Majorana fermions, where 
TQC by 'braiding' the vortices is possible. Several works 
have addressed properties of inter- vortex Majorana tun- 
neling both in the case of a single vortex pair[9j [10] as 
well as in vortex lattices [TTUT3] . In particular, for infinite 
vortex lattices with finite magnetic penetration depth, A, 
methods exist for accurately computing the quasiparti- 
cle spectrum [TT| [T2] . In this paper we shall consider the 
opposite limit where A is much larger than the lattice 
size. This condition can be satisfied for 2D SC films and 
is also the limit where TQC with Majorana fermions is 
best understood. 

In this work we analyze the problem of Majorana 
fermions tunneling within vortex lattices in chiral p-wave 
superconductors (or its equivalents [6 J and discover that 



the order parameter phase plays a crucial role, leading 
to surprising and promising features. We find that the 
purely imaginary tunneling amplitude tAB between the 
Majorana fermions at two vortices A and B depends on 
the order parameter phases a and /3 at the two vortices 
according to tAB oc is'm[(a — /3)/2]. This dependence is 
a feature of Majorana fermions since they require tAB — 
—tsA- The amplitude also decays exponentially and os- 
cillates with the distance between the vortices, in agree- 
ment with previous work with two vortices [9( [TQl HI] . 
Numerically computing nearest neighbor hopping am- 
plitudes in a large triangular/square lattice of vortices 
shows that in the central regions of these lattices and in 
particular, covering almost the entire lattice when it has a 
circular boundary, we find a superlattice emerging due to 
amplitudes vanishing on certain bonds. The Hamiltonian 
there is found to be translation-invariant with a super- 
lattice lattice constant that is twice that of the original 
vortex lattice, up to appropriate Z2 gauge transforma- 
tions. In both the square and triangular lattice cases we 
find localizable flat bands of Majorana fermions at zero 
energy, i.e, with no energy shift as compared with their 
energy in an isolated vortex. This macroscopic degener- 
acy of the ground state should dominate thermodynamic 
and transport characteristics of the Majorana fermions 
in these lattices and will also result in stronger quantum 
coherence than has been believed till now. We also find 
bands with non-zero energy which are gapped in the case 
of the triangular lattice and gapless for the square lattice. 

Inter-vortex tunneling amplitude: We consider a 
network of vortices centered at locations {r^} with the 
same vorticity L Following past research |1, we shall ne- 
glect the effects of the magnetic vector potential and as- 
sume that A is large compared to the lattice size. The 
chiral p-wave superconducting(SC) order parameter is 
given by A(r)(d x + id y ), where A(r) = f(r)e^ r \ The 
SC amplitude f(r) is real, positive and equal to Aq at 
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FIG. 1. (Color online) Triangular and square vortex lat- 
tices with circular boundaries showing the magnitudes and 
signs of the purely imaginary nearest neighbor hopping am- 
plitudes between Majorana states, encoded by the opacity 
and color (green/red=+/— ) of the directed bonds respec- 
tively. The translation-invariant Zi flux per plaquette is neg- 
ative/positive for the triangular/square lattice as there are 
an odd/even number of negative (red) directed bonds around 
each loop. Inset are zoomed sections showing the numbered 
four-site bases of the super lattices. 



all places other than within a distance of ~ £ (coherence 
length) from the vortex cores, where it drops to zero with 
a functional dependence fo(\r — Vk\) near the k th core. 
<j)(r) is the phase of the SC order parameter. We now 
isolate two vortices A and B, with vortex B located at 
a polar coordinate (7£, ©) with respect to vortex A. The 



polar coordinates of a generic point is denoted by (r, 0) 
or (p, d) with respect to A or B respectively. The (real) 
phase <p(r) takes the value a, f3 immediately to the 'right' 
of the vortices A and B respectively, so that near A and 
B the vortex phase (j)(r) varies as a + £0 and j3 + M 
respectively. 

The unperturbed zero energy Majorana wavefunctions 
near each of the two vortices, that exist only if the vor- 
ticity £ is odd, are given by (for a derivation see [15] for 
the £ = 1 case) 



1p A OC 



l/j B OC 



e ia/2 e i(l+l)6/2 
_ e -ia/2 e -i(l+l)6/2 

e i(3/2 e i(£+l)$/2 
_ e -i(3/2 e -i(£+l)$/2 






where 



$t(r) = J ( £ + i )/2 (^r)exp ( -—- / f (r')dr' 



(la) 
(lb) 

(2) 



and k, = kp \A — Aq/(^f&f) 2 - We shall assume that the 
lowest energy states in the presence of tunneling between 
the two vortices is spanned by these two wavefunctions. 
Introducing the notation Oa = {A, (d x + id y )} /2 and H 
the single-particle part of the effective Hamiltonian, the 
amplitude for tunneling from vortex A to B is given by: 
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i ff d 2 x Im [ e -i(-+P)/2^ p)e -i(i+i^/2 

x(0A-0 AA )e-^ +1 W 2 $(r)], (3) 

where the proportionality factors are real numbers and 
will be in the following equations. Aa introduced in the 
second step above is the order parameter when only vor- 
tex A is present and whose phase is given by a + £0, i.e, 
near vortex A, A ~ A^. Thus, the overlap integral be- 
comes significant only if we are far enough from A. Mean- 
while we see that the exponentially decaying factors in 
the functions <£ confine any meaningful contribution from 
the overlap integral to a strip of area connecting vortices 
A and B and having a width of the order of the coherence 
length £. Finally derivatives of A near vortex B will be 
larger by a factor of the order of (H/%) when compared 
to derivatives of ipA and A^ . Thus we estimate that 



t A B oc i (J d 2 xlm[e-^ a+ ^/ 2 ^(p)^(r) 

J J near B 

< e+1 ^ +e ^ 2 (d x + id y )A(r)]. 
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Near B, A(r) « f (p)e if3 e iM and so 

Thus, 
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J J near B 
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Since the above integrand transforms to its complex con- 
jugate when we consider points symmetric about the line 
joining the two vortices, the integral is real and so we 
arrive at a central result of this paper 



tAB oc z sin 



(7) 



The hopping amplitude is purely imaginary as it must 
be for Majorana states. The proportionality factor in- 
volves <E>(7£) and its derivatives which give rise to the ex- 
ponential decay and oscillatory behavior that have been 
obtained before [9 . We note here that there is an ambi- 
guity of 2tt in each of the phases a and /3 and this leads 
to an ambiguity in the sign of all amplitudes connecting 
a given vortex - this is however physically irrelevant due 
to the Z2 freedom in the sign of the Majorana operator. 
Such dependence on the SC phase that is anti-symmetric 
under the exchange of the vortices is essential because if 
the hopping amplitude were a function only of the scalar 
7£, then it has to be zero since it has to change sign when 
the two Majorana states are interchanged, i.e, we require 
tAB — —tBA- This physics seems to have been obscured 
in previous research on this topic. 

Tight-binding models: We consider nearest neighbor 
hopping only as by comparison all longer range hops are 
exponentially suppressed with distance. Numerical cal- 
culations in the following sections using the ansatz Eq. 
([7]) suggest simple translationally invariant superlattice 
models for the Majorana fermions when dealing with the 
triangular and square vortex lattices (similar calculations 
may be performed for other lattices). We consider odd 
vorticities £ since only in those cases can we have the Ma- 
jorana states and in those cases numerical computation 
yields a simple picture. The SC phase profile is calculated 
as(/>(r) = £^2 k SiTg[z(r)-z(r k )], where z(x,y) = x+iy$\. 
While the exact realization of the hopping amplitude pro- 
file depends sensitively on the shape of the lattice bound- 
ary, in the bulk central regions of all lattices and covering 
almost all of the area of a lattice with a circular boundary 
we find a superlattice network as shown in Figure [I] As 
can be seen there, many hopping amplitudes go to zero 
and the rest define superlattices with constant Z2 flux 
per plaquette. These should thus be gauge-equivalent to 
translation-invariant tight binding models for the Majo- 
rana fermions, as derived below. 
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(a)DOS for triangular lattice 
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FIG. 2. (Color online) Histograms of energy eigenstates in the 
lattices shown in Figure [I] Both show macroscopic degenera- 
cies at zero energy, with some broadening due to boundary 
effects. The density of states (DOS) at higher energies are fit 
well using energy bands calculated in the main text (Eqs ( 11 ) 
and (pt). 



The triangular lattice: The superlattice here is a 
triangular lattice with a basis of 4 vortices, as shown in 



Figure 1(a) The hopping amplitudes along the directed 
intra-basis bonds are the same while the inter-site am- 
plitudes have the opposite sign. The Hamiltonian for the 
bulk is thus (up to an overall scale factor): 



R a=l 



(8) 



The S a are marked on Figure 1(a) while 7^ is the Ma- 
jorana state operator for the a tn site at superlattice 
point R. We choose the usual convention normalizing 

the 7 operators S7^,7^/ \ = ^aa'^R,R' • The Hamilto- 
nian ([8| can be solved by using the traveling modes [16] 
b^ = —7= r ^2nJR e iq R , where J\f is the number of su- 
perlattice sites. The q vectors span the BZ of the su- 
perlattice and the modes at ±q are identified with each 
other: b% = (b°L q ) ] and {&£,&£} = S a(3 S q - p . The Hamil- 
tonian can now be rewritten as 

Kb = \T, J2 (^) + ^(9)6g, (9) 
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a = 0,/3 = 1,2,3, 

a = 1,2,3, (3 = 0, (10) 

otherwise. 



Finding the excitation spectrum of this model is re- 
duced to diagonalizing h(q) and retaining an indepen- 
dent half of the particle- hole conjugate modes with non- 
negative eigenvalues. In this case, diagonalization yields 
two bands of complex fermions with energies 



Ei(q) = 0, 



(11a) 



E2(q) = 8W 1 + 2 cos ^— - cos ^— - cos ^— -. (lib) 

The eigenfunctions corresponding to the zero energy flat 
band E\(q) also possess Chern number equal to zero, 
implying that we can form Wannier- localized Major ana 
states involving the bond center sites whose energies are 
(still) zero! As shown in Figure 2(a)| the DOS obtained 
from this model is a good approximation to the numeri- 



cally computed DOS of the model in Figure [1(a) 
The square lattice: The superlattice for this problem 
is a square lattice with a basis of 4 vortices, as shown 
in Figure 1(b) The hopping amplitudes along directed 
bonds can be chosen to be the same. The Hamiltonian 
for the bulk is thus (up to an overall scale factor): 



^^EE^K-^hs + h-c- 



(12) 



R a=l 



The basis sites labeled '3' are disconnected from the rest 
of the lattice and so they form a flat band of site-localized 
zero energy Majorana states! Using techniques intro- 
duced in the previous section, we find that the other 
three sites lead to the formation of another zero energy 
flat band of Majorana fermions as well as a gapless band 
of complex fermions spanning the entire BZ and having 
the dispersion 



E 3 (q) 



sm 2 (q x /2)+sm 2 (q y /2). 



(13) 



As in the case of the triangular lattice, the zero energy 
band wavefunctions have Chern number zero and so can 
be used to form zero energy Wannier- localized Majorana 
states. As shown in Figure [2(b)] the DOS obtained from 
this model is a good approximation to the numerically 



computed DOS of the lattice in Figure 1(b) 
Discussion and conclusions: We have discovered that 
the global superconducting phase profile plays an im- 
portant role in the tunneling amplitudes of Majorana 
fermions tunneling between vortices. This is necessary 
since the hopping amplitude tAB needs to satisfy the re- 
quirement tAB — —tBA for Majorana fermions, which is 
impossible if tAB is a nonzero and a function only of the 
inter-vortex distance. This conclusion was obscured in 



previous research involving only two vortices, for which 
the relative phase \a — /3\ = £tt is untunable. We estimate 
that the SC phase dependence of the tunneling amplitude 
manifests itself in the sinusoidal law expressed in equa- 
tion ([7]). This has surprising consequences for triangular 
and square vortex lattices where we find emergent su- 
perlattice structures hosting zero energy localizable Ma- 
jorana fermion states that number a finite fraction of the 
original number of states. In the square lattice these also 
include single-site localized states on a quarter of all lat- 
tice sites. These zero energy states imply the existence of 
temporally coherent states even when the vortices are rel- 
atively close together which bodes well for the feasibility 
of TQC in such setups. 

We emphasize here that the design and implementation 
of such vortex lattices is feasible [17] and not just an es- 
oteric possibility. Along the line of systems that have 
already been fabricated [8 , such vortex lattices may be 
created by depositing a layer of s-wave superconducting 
metal like Nb or Al on a strong topological insulator(STI) 
such as Bi2Se3. The deposited layer is continuous except 
for the presence of holes laid out in the pattern of the 
required vortex lattice. Cooling this arrangement below 
a suitable temperature will induce superconductivity in 
the STI surface state that has been predicted to host 
Majorana states at vortices [6 . Turning on a magnetic 
field of sufficient strength will cause vortices to appear 
at the holes since the induced superconductor is weak- 
est at those points, providing a strong pinning potential. 
Designing protocols for using these zero energy states in 
the vortex lattice for TQC is an open problem. 
An important ingredient in obtaining the phase- 
dependence of the tunneling amplitude is the choice of 
trial wavefunctions (Eq. (fTl) spanning the low energy 
subspace, which is the same as in other previous ap- 
proaches (e.g [9]). These wavefunctions may be improved 
to take into account the mismatch between the phase of 
A(r) and the uniform angular twisting of the SC order 
parameter phase, in the spirit of [TTJ[T2] but modified for 
our scenario where the magnetic length is infinite. 
Finally, while considering disorder in the tunneling am- 
plitudes that arise due to randomness in the inter- vortex 
distances, our work suggests that past calculations [T8H22] 
be extended to take into account the long range effects 
of lattice geometry on the effective anyonic Hamiltonian, 
mediated by the order parameter phase profile. For e.g, 
in the case mentioned here the bonds with zero strength 
have negligible noise and do not fit into the scheme of ear- 
lier analyses. Also, we consider here lattices with circular 
boundary which are very well-approximated by transla- 
tion invariant Majorana tight-binding models. The case 
of a generic boundary - a candidate for future studies - 
is different as there is an annular region with apparently 
disordered bonds but which still leave a finite fraction 
of bulk Majorana states at zero energy. In closing, we 
note that understanding the effects of including the gauge 



field [TTJ [12] as well as the effect of SC phase profiles that 
are not Abrikosov lattice-like are also challenging avenues 
of future research. 
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